Load the necessary libraries
library(car) #for regression diagnostics
library(broom) #for tidy output
library(broom.mixed)
library(ggfortify) #for model diagnostics
library(sjPlot) #for outputs
library(knitr) #for kable
library(effects) #for partial effects plots
library(ggeffects) #for effects plots in ggplot
library(emmeans) #for estimating marginal means
library(MASS) #for glm.nb
library(MuMIn) #for AICc
library(tidyverse) #for data wrangling
library(DHARMa) #for residuals and diagnostics
library(nlme) #for lme
library(lme4) #for glmer
library(glmmTMB) #for glmmTMB
library(performance) #for diagnostic plots
library(see) #for diagnostic plots
To investigate synergistic coral defence by mutualist crustaceans, Mckeon et al. (2012) conducted an aquaria experiment in which colonies of a coral species were placed in a tank along with a preditory seastar and one of four symbiont combinations:
The experiments were conducted in a large octagonal flow-through seawater tank that was partitioned into eight sections, which thereby permitted two of each of the four symbiont combinations to be observed concurrently. The tank was left overnight and in the morning, the presence of feeding scars on each coral colony was scored as evidence of predation. The experiments were repeated ten times, each time with fresh coral colonies, seastars and symbiont.
The ten experimental times represent blocks (random effects) within which the symbiont type (fixed effect) are nested.
mckeon = read_csv('../public/data/mckeon.csv', trim_ws=TRUE)
##
## ── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────
## cols(
## BLOCK = col_double(),
## PREDATION = col_double(),
## SYMBIONT = col_character()
## )
glimpse(mckeon)
## Rows: 80
## Columns: 3
## $ BLOCK <dbl> 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10…
## $ PREDATION <dbl> 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, …
## $ SYMBIONT <chr> "none", "none", "none", "none", "none", "none", "none", "non…
Since the response here is the presence or absence of predation (feeding scars), a binomial distribution is appropriate.
We need to make sure that the categorical variable and the random effect are defined as factors. When doing so, it might be valuable to rearrange the order of the fixed effect (SYMBIONT) such that the none group is considered the first group. This way, the other levels will all naturally be compared to this level (hence it will be treated as a reference of control group).
mckeon = mckeon %>%
mutate(BLOCK=factor(BLOCK),
SYMBIONT=factor(SYMBIONT, levels=c('none','crabs','shrimp','both')))
Model formula: \[ y_i \sim{} \mathcal{N}(n, p_i)\\ ln\left(\frac{p_i}{1-p_1}\right) =\boldsymbol{\beta} \bf{X_i} + \boldsymbol{\gamma} \bf{Z_i} \]
where \(\boldsymbol{\beta}\) and \(\boldsymbol{\gamma}\) are vectors of the fixed and random effects parameters respectively and \(\bf{X}\) is the model matrix representing the overall intercept and effects of symbionts on the probability of the colony experiencing predation. \(\bf{Z}\) represents a cell means model matrix for the random intercepts associated with individual coral colonies.
ggplot(mckeon, aes(y=PREDATION, x=SYMBIONT)) +
geom_point(position=position_jitter(width=0.2, height=0))+
facet_wrap(~BLOCK)
The data have been setup as a single factor (mixed effects) design. Alternatively, we could have coded this up as a factorial (mixed effects) such that there was a dummy variable for Crabs, and a dummy variable for Shrimps. We could then fit a model that had the main effects of Crabs and Shrimp as well as there interaction.
Since models with categorical variables essentially dummy code up the variable anyway, fitting a single factor model will result in similar parameter estimates to a factorial design with the exception of the interaction (or Crabs and Shrimp group). Nevertheless, having fit the model as a single factor model, we can then perform a specific contrast that is the same as the interaction effect from a factorial design.
Lets start with the random intercept model
mckeon.glmer1 = glmer(PREDATION ~ SYMBIONT+(1|BLOCK), data=mckeon,
family=binomial(link='logit'))
And now attempt the random intercept/slope model. By default, glmer uses the bobyqa optimizer. This optimizer uses iterative quadratic approximation of twice differentiable functions. When functions are not twice differentiable, it may perform poorly.
We could specify the new random intercept/slope model out in full or we could simply use the update() function to make a modification (in this case, we can remove the random intercept and add the random intercept/slope).
mckeon.glmer2 <- update(mckeon.glmer1, ~ . - (1|BLOCK) + (SYMBIONT|BLOCK))
## Warning in (function (fn, par, lower = rep.int(-Inf, n), upper = rep.int(Inf, :
## failure to converge in 10000 evaluations
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge: degenerate Hessian with 1 negative eigenvalues
Along with attempting to fit the model, a glmer() performs a set of model checks to help diagnose whether the optimization process has successfully converged as well as passed other criterion.
Evidently, the above (random intercepts/slope) model failed to converge after 10000 evaluations of the likelihood function. That is, it failed to reach a consensus on the parameter estimates. Note, this is a warning not an error. It does not necessarily mean that the model fit is not correct. Nevertheless, without further exploration of what might be causing the lack of convergence, such a model should not be used as the paramter estimates (and/or their standard errors) may be unreliable.
When receiving a convergence warning, the following potential remedies are recommended:
alternative optimizers (such as Nelder-Mead, nlminbwrap or L-BFGS-B)
increase the number of iteratations in an attempt to provide more scope for consensus
abandon the model
use allFit() to attempt to fit will all the available optimizers. Different optimizers perform differently under different conditions. The default optimizer (bobyqa) has been selected by the lme4 team as they consider it to be the best all-round choice and a good compromise between robustness and efficiency. However, other optimizers might work better for any given circumstance. Other common optimizers include:
If all of the optimizers yield very similar parameter estimates, then the convergence warnings can be considered as false positives and as such, ignored (unless they all suffer the same structural issue due to a poorly speficied or overfit model).
adjust the stopping tolerance for the optimizer (e.g. optCtrl=list(maxfun=2e5))
center (and possibly scale) continuous predictors
consider a more expensive Hessian calculation. Hessian compuations are used for convergence checking as well as estimating the standard errors of fixed effects parameters in GLMM. By default, the method used provides very rapid estimates. However, they can be unreliable.
restart the fit with starting values taken from the current parameter estimates
Let start with the first of these options - often it is the best initial approach.
mckeon.allFit <- allFit(mckeon.glmer2)
## Loading required namespace: dfoptim
## Loading required namespace: optimx
## bobyqa :
## boundary (singular) fit: see ?isSingular
## [OK]
## Nelder_Mead :
## boundary (singular) fit: see ?isSingular
## [OK]
## nlminbwrap :
## boundary (singular) fit: see ?isSingular
## [OK]
## nmkbw : [OK]
## optimx.L-BFGS-B :
## Warning in optwrap(optimizer, devfun, start, rho$lower, control = control, :
## convergence code 52 from optimx
## Warning in optimx.check(par, optcfg$ufn, optcfg$ugr, optcfg$uhess, lower, : Parameters or bounds appear to have different scalings.
## This can cause poor performance in optimization.
## It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.
## boundary (singular) fit: see ?isSingular
## [OK]
## nloptwrap.NLOPT_LN_NELDERMEAD :
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge: degenerate Hessian with 1 negative eigenvalues
## [OK]
## nloptwrap.NLOPT_LN_BOBYQA :
## boundary (singular) fit: see ?isSingular
## [OK]
## Check which of the models are considered valid (OK)
is.OK <- sapply(mckeon.allFit, is, "merMod")
is.OK
## bobyqa Nelder_Mead
## TRUE TRUE
## nlminbwrap nmkbw
## TRUE FALSE
## optimx.L-BFGS-B nloptwrap.NLOPT_LN_NELDERMEAD
## TRUE TRUE
## nloptwrap.NLOPT_LN_BOBYQA
## TRUE
Conclusions:
We will check the convergence characteristics of only the valid models
diff_optims.OK <- mckeon.allFit[is.OK]
lapply(diff_optims.OK,function(x) x@optinfo$conv$lme4$messages)
## $bobyqa
## [1] "boundary (singular) fit: see ?isSingular"
##
## $Nelder_Mead
## [1] "boundary (singular) fit: see ?isSingular"
##
## $nlminbwrap
## [1] "boundary (singular) fit: see ?isSingular"
##
## $`optimx.L-BFGS-B`
## [1] "boundary (singular) fit: see ?isSingular"
##
## $nloptwrap.NLOPT_LN_NELDERMEAD
## [1] "unable to evaluate scaled gradient"
## [2] "Model failed to converge: degenerate Hessian with 1 negative eigenvalues"
##
## $nloptwrap.NLOPT_LN_BOBYQA
## [1] "boundary (singular) fit: see ?isSingular"
Conclusions:
Singularity refers to presence of one or more parameters whose estimated values are on the boundary of the possible parameter space (range of possible values). Typically, this is either that there are either variances that are estimated to be 0, or correlations that are estimated to be either -1 or 1.
Singular fits usually suggest overfitting and are more likely to contain numerical problems and lead to issues concerning testing hypotheses on boundaries.
Situations that are likely to contribute to singularity include:
f|g - particularly where f is a categorical variable with a relatively large number of levels)Advice on dealing with singlularity issues varies widely from simplifying the random effects structure to fitting the model in a Bayesian framework in which singularity can be avoided via appropriate priors.
For the current example, we might consider dropping the random intercept/slope.
Before excepting that we must drop the random intercept/slope, it would be worth quickly exploring the other options listed above.
Lets see if extending the number of evaluations will result in convergence.
mckeon.glmer2 <- update(mckeon.glmer1, ~ . - (1|BLOCK) + (SYMBIONT|BLOCK),
control=glmerControl(optimizer='bobyqa',
optCtrl=list(maxfun=2e5)))
## boundary (singular) fit: see ?isSingular
Conclusions:
In the current example, the fixed effect (SYMBIONT) is a categorical variable. Consequently, it is effectively already scaled, so the centering (and/or scaling) option isn’t a solution here.
Now lets consider the more expensive Hessian computation.
pars <- getME(mckeon.glmer2, c("theta","fixef"))
devfun <- update(mckeon.glmer2, devFunOnly=TRUE)
library(numDeriv)
cat("hess:\n"); print(hess <- hessian(devfun, unlist(pars)))
## hess:
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.4213975722 -2.026547e-04 -1.644069e-01 1.526993e-01 -2.826724e-04
## [2,] -0.0002026547 1.473090e-05 -1.103229e-04 2.360222e-07 5.889714e-05
## [3,] -0.1644068830 -1.103229e-04 5.610199e-01 -4.732408e-01 -3.033088e-04
## [4,] 0.1526992889 2.360222e-07 -4.732408e-01 4.070909e-01 -7.901643e-05
## [5,] -0.0002826724 5.889714e-05 -3.033088e-04 -7.901643e-05 2.743109e-04
## [6,] 0.0595859600 -3.033426e-04 -1.002367e-01 1.109353e-01 -1.446770e-03
## [7,] -0.0128592645 -7.902433e-05 1.109353e-01 -8.059782e-02 -3.480331e-04
## [8,] -0.0001698081 -9.029730e-09 -1.602309e-07 -3.028307e-07 -3.996483e-09
## [9,] -0.0010049319 -2.569394e-08 -6.284820e-07 -4.679696e-07 -1.400572e-08
## [10,] -0.0009648772 -5.504486e-08 -6.568830e-07 -3.542815e-07 -1.528416e-08
## [11,] 0.0017816129 8.246683e-04 -1.147059e-02 -2.284188e-03 3.813384e-03
## [12,] 0.0005300320 -9.729221e-05 5.407133e-04 1.195055e-04 -4.498222e-04
## [13,] -0.0353598093 7.565127e-04 -3.488098e-02 -8.343249e-03 3.498610e-03
## [14,] 0.0366325791 1.654481e-04 2.286968e-02 5.939556e-03 7.655850e-04
## [,6] [,7] [,8] [,9] [,10]
## [1,] 5.958596e-02 -1.285926e-02 -1.698081e-04 -1.004932e-03 -9.648772e-04
## [2,] -3.033426e-04 -7.902433e-05 -9.029730e-09 -2.569394e-08 -5.504486e-08
## [3,] -1.002367e-01 1.109353e-01 -1.602309e-07 -6.284820e-07 -6.568830e-07
## [4,] 1.109353e-01 -8.059782e-02 -3.028307e-07 -4.679696e-07 -3.542815e-07
## [5,] -1.446770e-03 -3.480331e-04 -3.996483e-09 -1.400572e-08 -1.528416e-08
## [6,] 1.192637e-01 1.566713e-02 -1.551708e-07 -1.032991e-07 -1.750863e-07
## [7,] 1.566713e-02 5.188491e-02 9.812057e-08 -2.048719e-07 -1.002466e-07
## [8,] -1.551708e-07 9.812057e-08 5.836355e-01 -4.984751e-01 -1.088363e-03
## [9,] -1.032991e-07 -2.048719e-07 -4.984751e-01 4.260798e-01 -9.365976e-04
## [10,] -1.750863e-07 -1.002466e-07 -1.088363e-03 -9.365976e-04 4.259329e-01
## [11,] -5.303501e-02 -1.056178e-02 -9.029329e-08 -4.217722e-07 -4.492473e-07
## [12,] 2.500050e-03 5.533505e-04 -1.728246e-08 -2.747050e-08 -1.480038e-08
## [13,] -1.612714e-01 -3.857528e-02 -2.337290e-07 -6.183609e-07 -5.379129e-07
## [14,] 1.057451e-01 2.746087e-02 -7.623707e-09 -8.387474e-08 -1.079613e-07
## [,11] [,12] [,13] [,14]
## [1,] 1.781613e-03 5.300320e-04 -3.535981e-02 3.663258e-02
## [2,] 8.246683e-04 -9.729221e-05 7.565127e-04 1.654481e-04
## [3,] -1.147059e-02 5.407133e-04 -3.488098e-02 2.286968e-02
## [4,] -2.284188e-03 1.195055e-04 -8.343249e-03 5.939556e-03
## [5,] 3.813384e-03 -4.498222e-04 3.498610e-03 7.655850e-04
## [6,] -5.303501e-02 2.500050e-03 -1.612714e-01 1.057451e-01
## [7,] -1.056178e-02 5.533505e-04 -3.857528e-02 2.746087e-02
## [8,] -9.029329e-08 -1.728246e-08 -2.337290e-07 -7.623707e-09
## [9,] -4.217722e-07 -2.747050e-08 -6.183609e-07 -8.387474e-08
## [10,] -4.492473e-07 -1.480038e-08 -5.379129e-07 -1.079613e-07
## [11,] 3.960248e+00 -6.832716e-03 3.526121e-01 1.446899e-02
## [12,] -6.832716e-03 7.879640e-04 -6.203661e-03 -1.417017e-03
## [13,] 3.526121e-01 -6.203661e-03 1.164430e+00 -8.056139e-01
## [14,] 1.446899e-02 -1.417017e-03 -8.056139e-01 8.214999e-01
cat("grad:\n"); print(grad <- grad(devfun, unlist(pars)))
## grad:
## [1] 3.633557e-03 -7.390718e-09 3.913262e-08 -5.231554e-08 -4.078292e-08
## [6] 4.978622e-09 1.356960e-09 -6.296769e-08 5.468531e-08 1.552051e-08
## [11] 1.583982e-07 -2.453373e-08 -2.122501e-08 -5.751713e-08
cat("scaled gradient:\n")
## scaled gradient:
#print(scgrad <- solve(chol(hess), grad))
## compare with internal calculations:
mckeon.glmer2@optinfo$derivs
## $gradient
## [1] 3.633557e-03 -7.105427e-09 4.288125e-08 -5.375256e-08 -4.074963e-08
## [6] 4.476419e-09 1.350031e-09 -6.313172e-08 5.453415e-08 1.463718e-08
## [11] 1.515588e-07 -2.437162e-08 -2.149392e-08 -5.524470e-08
##
## $Hessian
## [,1] [,2] [,3] [,4] [,5]
## [1,] 4.206285e-01 -2.052784e-04 -1.644053e-01 1.526983e-01 -2.810955e-04
## [2,] -2.052784e-04 2.288818e-05 -1.113415e-04 7.152557e-07 5.769730e-05
## [3,] -1.644053e-01 -1.113415e-04 5.610199e-01 -4.732451e-01 -3.056526e-04
## [4,] 1.526983e-01 7.152557e-07 -4.732451e-01 4.070864e-01 -7.963181e-05
## [5,] -2.810955e-04 5.769730e-05 -3.056526e-04 -7.963181e-05 2.880096e-04
## [6,] 5.958509e-02 -3.020763e-04 -1.002369e-01 1.109383e-01 -1.446962e-03
## [7,] -1.285934e-02 -7.915497e-05 1.109328e-01 -8.060026e-02 -3.476143e-04
## [8,] 2.384186e-07 1.192093e-06 2.384186e-07 2.384186e-07 -7.152557e-07
## [9,] 1.668930e-06 -1.907349e-06 -4.768372e-07 -2.384186e-07 9.536743e-07
## [10,] 1.430511e-06 2.384186e-07 1.907349e-06 1.907349e-06 9.536743e-07
## [11,] 1.781225e-03 8.265972e-04 -1.146817e-02 -2.286196e-03 3.813505e-03
## [12,] 5.316734e-04 -9.822845e-05 5.424023e-04 1.187325e-04 -4.477501e-04
## [13,] -3.535891e-02 7.569790e-04 -3.487992e-02 -8.344889e-03 3.500223e-03
## [14,] 3.663301e-02 1.678467e-04 2.287388e-02 5.942345e-03 7.662773e-04
## [,6] [,7] [,8] [,9] [,10]
## [1,] 5.958509e-02 -1.285934e-02 2.384186e-07 1.668930e-06 1.430511e-06
## [2,] -3.020763e-04 -7.915497e-05 1.192093e-06 -1.907349e-06 2.384186e-07
## [3,] -1.002369e-01 1.109328e-01 2.384186e-07 -4.768372e-07 1.907349e-06
## [4,] 1.109383e-01 -8.060026e-02 2.384186e-07 -2.384186e-07 1.907349e-06
## [5,] -1.446962e-03 -3.476143e-04 -7.152557e-07 9.536743e-07 9.536743e-07
## [6,] 1.192703e-01 1.566839e-02 9.536743e-07 0.000000e+00 -7.152557e-07
## [7,] 1.566839e-02 5.189705e-02 1.192093e-06 -1.907349e-06 2.384186e-07
## [8,] 9.536743e-07 1.192093e-06 5.827065e-01 -4.972365e-01 -2.384186e-07
## [9,] 0.000000e+00 -1.907349e-06 -4.972365e-01 4.245358e-01 -1.668930e-06
## [10,] -7.152557e-07 2.384186e-07 -2.384186e-07 -1.668930e-06 4.245272e-01
## [11,] -5.303502e-02 -1.056194e-02 -4.768372e-07 -4.768372e-07 -4.768372e-07
## [12,] 2.501011e-03 5.519390e-04 -2.384186e-06 2.622604e-06 7.152557e-07
## [13,] -1.612737e-01 -3.857756e-02 -2.384186e-07 2.384186e-06 1.192093e-06
## [14,] 1.057370e-01 2.746010e-02 -1.907349e-06 2.145767e-06 -4.768372e-07
## [,11] [,12] [,13] [,14]
## [1,] 1.781225e-03 5.316734e-04 -3.535891e-02 3.663301e-02
## [2,] 8.265972e-04 -9.822845e-05 7.569790e-04 1.678467e-04
## [3,] -1.146817e-02 5.424023e-04 -3.487992e-02 2.287388e-02
## [4,] -2.286196e-03 1.187325e-04 -8.344889e-03 5.942345e-03
## [5,] 3.813505e-03 -4.477501e-04 3.500223e-03 7.662773e-04
## [6,] -5.303502e-02 2.501011e-03 -1.612737e-01 1.057370e-01
## [7,] -1.056194e-02 5.519390e-04 -3.857756e-02 2.746010e-02
## [8,] -4.768372e-07 -2.384186e-06 -2.384186e-07 -1.907349e-06
## [9,] -4.768372e-07 2.622604e-06 2.384186e-06 2.145767e-06
## [10,] -4.768372e-07 7.152557e-07 1.192093e-06 -4.768372e-07
## [11,] 3.960271e+00 -6.834507e-03 3.526089e-01 1.446557e-02
## [12,] -6.834507e-03 7.982254e-04 -6.203413e-03 -1.418829e-03
## [13,] 3.526089e-01 -6.203413e-03 1.164444e+00 -8.056104e-01
## [14,] 1.446557e-02 -1.418829e-03 -8.056104e-01 8.215151e-01
Conclusions:
Finally, lets consider re-fitting the model using the current parameter estimates as starting values.
#strict_tol <- glmerControl(optCtrl=list(xtol_abs=1e-8, ftol_abs=1e-8))
pars <- getME(mckeon.glmer2, c("theta","fixef"))
mckeon.glmer3 <- update(mckeon.glmer2, start=pars, control=glmerControl(optimizer='bobyqa',
optCtrl=list(maxfun=2e5)))
## boundary (singular) fit: see ?isSingular
Conclusions:
Lets start with the random intercept model
mckeon.glmmTMB1 = glmmTMB(PREDATION ~ SYMBIONT+(1|BLOCK), data=mckeon,
family=binomial(link='logit'), REML=TRUE)
And now attempt the random intercept/slope model.
mckeon.glmmTMB2 <- update(mckeon.glmmTMB1, ~ . - (1|BLOCK) + (SYMBIONT|BLOCK))
## Error in (function (start, objective, gradient = NULL, hessian = NULL, : NA/NaN gradient evaluation
## Timing stopped at: 0.099 0.001 0.099
mckeon.glmmTMB2 <- update(mckeon.glmmTMB1, ~ . - (1|BLOCK) + (SYMBIONT|BLOCK),
control=glmmTMBControl(optimizer=optim,
optArgs = list(method='SANN')))
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')
AICc(mckeon.glmmTMB1, mckeon.glmmTMB2)
IN addition to numerous warnings, there is an error about NA/NaN gradient evaluation. This is a somewhat generic error returned by the nlminb optimizer. It occurs when there is an error in either the objective function or gradient function calculations.
Ultimately, this might suggest a poorly specified or overfit model. We will drop the random intercept/slope.
plot_grid(plot_model(mckeon.glmer1, type='diag'))
Conclusions:
performance::check_model(mckeon.glmer1)
Conclusions:
mckeon.resid = simulateResiduals(mckeon.glmer1, plot=TRUE)
testZeroInflation(mckeon.resid)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0428, p-value = 0.936
## alternative hypothesis: two.sided
Conclusions:
plot_grid(plot_model(mckeon.glmmTMB1, type='diag'))
Conclusions:
performance::check_model(mckeon.glmmTMB1)
Conclusions:
mckeon.resid = simulateResiduals(mckeon.glmmTMB1, plot=TRUE)
testZeroInflation(mckeon.resid)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 0.96105, p-value = 0.968
## alternative hypothesis: two.sided
testDispersion(mckeon.resid)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## ratioObsSim = 1.0036, p-value = 0.928
## alternative hypothesis: two.sided
testZeroInflation(mckeon.resid)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 0.96105, p-value = 0.968
## alternative hypothesis: two.sided
#testTemporalAutocorrelation(mckeon.resid)
Conclusions:
plot_model(mckeon.glmer1, type='eff')
## $SYMBIONT
plot(allEffects(mckeon.glmer1))
ggpredict(mckeon.glmer1) %>% plot
## $SYMBIONT
ggemmeans(mckeon.glmer1, ~SYMBIONT) %>% plot
plot_model(mckeon.glmmTMB1, type='eff')
## $SYMBIONT
plot(allEffects(mckeon.glmmTMB1))
ggpredict(mckeon.glmmTMB1) %>% plot
## $SYMBIONT
ggemmeans(mckeon.glmmTMB1, ~SYMBIONT) %>% plot
summary(mckeon.glmer1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: PREDATION ~ SYMBIONT + (1 | BLOCK)
## Data: mckeon
##
## AIC BIC logLik deviance df.resid
## 70.7 82.6 -30.4 60.7 75
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -23.2076 -0.1730 0.0976 0.2980 0.8944
##
## Random effects:
## Groups Name Variance Std.Dev.
## BLOCK (Intercept) 11.81 3.437
## Number of obs: 80, groups: BLOCK, 10
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.096 1.811 2.814 0.00490 **
## SYMBIONTcrabs -3.842 1.465 -2.623 0.00871 **
## SYMBIONTshrimp -4.431 1.551 -2.856 0.00429 **
## SYMBIONTboth -5.599 1.724 -3.247 0.00116 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) SYMBIONTc SYMBIONTs
## SYMBIONTcrb -0.625
## SYMBIONTshr -0.652 0.736
## SYMBIONTbth -0.687 0.732 0.767
Conclusions:
tidy(mckeon.glmer1, effect='fixed', conf.int=TRUE, infer=c(TRUE, TRUE))
tidy(mckeon.glmer1, effect='fixed', conf.int=TRUE, infer=c(TRUE, TRUE), exponentiate=TRUE)
tidy(mckeon.glmer1, effect='fixed', conf.int=TRUE, infer=c(TRUE, TRUE), exponentiate=TRUE) %>% kable
| effect | term | estimate | std.error | statistic | p.value | conf.low | conf.high |
|---|---|---|---|---|---|---|---|
| fixed | (Intercept) | 163.3256181 | 295.7986843 | 2.813623 | 0.0048987 | 4.6929398 | 5684.1252579 |
| fixed | SYMBIONTcrabs | 0.0214459 | 0.0314110 | -2.623280 | 0.0087088 | 0.0012151 | 0.3785015 |
| fixed | SYMBIONTshrimp | 0.0119023 | 0.0184657 | -2.856063 | 0.0042893 | 0.0005689 | 0.2490134 |
| fixed | SYMBIONTboth | 0.0037016 | 0.0063822 | -3.247315 | 0.0011650 | 0.0001261 | 0.1086478 |
Conclusions:
# warning this is only appropriate for html output
sjPlot::tab_model(mckeon.glmer1, show.se=TRUE, show.aic=TRUE)
| PREDATION | ||||
|---|---|---|---|---|
| Predictors | Odds Ratios | std. Error | CI | p |
| (Intercept) | 163.33 | 295.80 | 4.69 – 5684.13 | 0.005 |
| SYMBIONT [crabs] | 0.02 | 0.03 | 0.00 – 0.38 | 0.009 |
| SYMBIONT [shrimp] | 0.01 | 0.02 | 0.00 – 0.25 | 0.004 |
| SYMBIONT [both] | 0.00 | 0.01 | 0.00 – 0.11 | 0.001 |
| Random Effects | ||||
| σ2 | 3.29 | |||
| τ00 BLOCK | 11.81 | |||
| ICC | 0.78 | |||
| N BLOCK | 10 | |||
| Observations | 80 | |||
| Marginal R2 / Conditional R2 | 0.228 / 0.832 | |||
| AIC | 70.706 | |||
summary(mckeon.glmmTMB1)
## Family: binomial ( logit )
## Formula: PREDATION ~ SYMBIONT + (1 | BLOCK)
## Data: mckeon
##
## AIC BIC logLik deviance df.resid
## 62.9 74.8 -26.4 52.9 79
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## BLOCK (Intercept) 17.83 4.223
## Number of obs: 80, groups: BLOCK, 10
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.420 1.782 2.481 0.013109 *
## SYMBIONTcrabs -3.317 1.313 -2.526 0.011534 *
## SYMBIONTshrimp -3.956 1.416 -2.793 0.005218 **
## SYMBIONTboth -5.152 1.565 -3.291 0.000997 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Conclusions:
tidy(mckeon.glmmTMB1, effect='fixed', conf.int=TRUE, infer=c(TRUE, TRUE))
## Warning in confint.glmmTMB(x, method = tolower(conf.method), level =
## conf.level, : extra arguments ignored: infer
tidy(mckeon.glmmTMB1, effect='fixed', conf.int=TRUE, infer=c(TRUE, TRUE), exponentiate=TRUE)
## Warning in confint.glmmTMB(x, method = tolower(conf.method), level =
## conf.level, : extra arguments ignored: infer
Conclusions:
# warning this is only appropriate for html output
sjPlot::tab_model(mckeon.glmmTMB1, show.se=TRUE, show.aic=TRUE)
| PREDATION | ||||
|---|---|---|---|---|
| Predictors | Odds Ratios | std. Error | CI | p |
| (Intercept) | 83.10 | 148.05 | 2.53 – 2730.03 | 0.013 |
| SYMBIONT [crabs] | 0.04 | 0.05 | 0.00 – 0.48 | 0.012 |
| SYMBIONT [shrimp] | 0.02 | 0.03 | 0.00 – 0.31 | 0.005 |
| SYMBIONT [both] | 0.01 | 0.01 | 0.00 – 0.12 | 0.001 |
| Random Effects | ||||
| σ2 | 3.29 | |||
| τ00 BLOCK | 17.83 | |||
| ICC | 0.84 | |||
| N BLOCK | 10 | |||
| Observations | 80 | |||
| Marginal R2 / Conditional R2 | 0.149 / 0.867 | |||
| AIC | 62.861 | |||
In addition to comparing each of the symbiont types against the control group of no symbionts, it might be interesting to investigate whether there are any differences between the predation protection provided by crabs and shrimp, as well as whether having both crabs and shrimp symbionts is different to only a single symbiont type.
These contrasts can be explored via specific contrasts.
| SYMBIONT | Crab vs Shrimp | One vs Both | None vs Symbiont |
|---|---|---|---|
| none | 0 | 0 | 1 |
| crab | 1 | 1/2 | -1/3 |
| shrimp | -1 | 1/2 | -1/3 |
| both | 0 | -1 | -1/3 |
cmat=cbind(
crab_vs_shrimp=c(0,1,-1,0),
one_vs_both=c(0,1/2,1/2,-1),
symbiont=c(1, -1/3, -1/3,-1/3)
)
round(crossprod(cmat),1)
## crab_vs_shrimp one_vs_both symbiont
## crab_vs_shrimp 2 0.0 0.0
## one_vs_both 0 1.5 0.0
## symbiont 0 0.0 1.3
## all contrasts orthogonal
emmeans(mckeon.glmer1, ~SYMBIONT, contr=list(SYMBIONT=cmat), type='link')
## $emmeans
## SYMBIONT emmean SE df asymp.LCL asymp.UCL
## none 5.096 1.81 Inf 1.55 8.65
## crabs 1.254 1.45 Inf -1.59 4.10
## shrimp 0.665 1.42 Inf -2.12 3.45
## both -0.503 1.40 Inf -3.25 2.24
##
## Results are given on the logit (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df z.ratio p.value
## SYMBIONT.crab_vs_shrimp 0.589 1.10 Inf 0.536 0.5921
## SYMBIONT.one_vs_both 1.462 1.02 Inf 1.429 0.1530
## SYMBIONT.symbiont 4.624 1.44 Inf 3.211 0.0013
##
## Results are given on the log odds ratio (not the response) scale.
emmeans(mckeon.glmer1, ~SYMBIONT, contr=list(SYMBIONT=cmat), type='response')
## $emmeans
## SYMBIONT prob SE df asymp.LCL asymp.UCL
## none 0.994 0.011 Inf 0.8243 1.000
## crabs 0.778 0.251 Inf 0.1690 0.984
## shrimp 0.660 0.319 Inf 0.1069 0.969
## both 0.377 0.329 Inf 0.0373 0.904
##
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale
##
## $contrasts
## contrast odds.ratio SE df z.ratio p.value
## SYMBIONT.crab_vs_shrimp 1.80 1.98 Inf 0.536 0.5921
## SYMBIONT.one_vs_both 4.32 4.42 Inf 1.429 0.1530
## SYMBIONT.symbiont 101.91 146.77 Inf 3.211 0.0013
##
## Tests are performed on the log odds ratio scale
Conclusions:
r.squaredGLMM(mckeon.glmer1)
## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help page.
## Warning: The null model is correct only if all variables used by the original
## model remain unchanged.
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0307141 (tol = 0.002, component 1)
## R2m R2c
## theoretical 0.2282003 0.8318487
## delta 0.2166278 0.7896639
## The delta method can be used with for all distributions and link
## functions, while lognormal approximation and trigamma function are
## limited to distributions with logarithmic link. Trigamma-estimate
## is recommended whenever available. Additionally, for binomial
## distributions, theoretical variances exist specific for each link
## function distribution.
performance::r2_nakagawa(mckeon.glmer1)
## # R2 for Mixed Models
##
## Conditional R2: 0.832
## Marginal R2: 0.228
cmat=cbind(
crab_vs_shrimp=c(0,1,-1,0),
one_vs_both=c(0,1/2,1/2,-1),
symbiont=c(1, -1/3, -1/3,-1/3)
)
round(crossprod(cmat),1)
## crab_vs_shrimp one_vs_both symbiont
## crab_vs_shrimp 2 0.0 0.0
## one_vs_both 0 1.5 0.0
## symbiont 0 0.0 1.3
## all contrasts orthogonal
emmeans(mckeon.glmmTMB1, ~SYMBIONT, contr=list(SYMBIONT=cmat), type='link')
## $emmeans
## SYMBIONT emmean SE df lower.CL upper.CL
## none 4.420 1.78 79 0.874 7.97
## crabs 1.103 1.61 79 -2.111 4.32
## shrimp 0.464 1.60 79 -2.719 3.65
## both -0.732 1.58 79 -3.871 2.41
##
## Results are given on the logit (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## SYMBIONT.crab_vs_shrimp 0.639 1.15 79 0.556 0.5800
## SYMBIONT.one_vs_both 1.516 1.04 79 1.454 0.1499
## SYMBIONT.symbiont 4.142 1.26 79 3.276 0.0016
##
## Results are given on the log odds ratio (not the response) scale.
emmeans(mckeon.glmmTMB1, ~SYMBIONT, contr=list(SYMBIONT=cmat), type='response')
## $emmeans
## SYMBIONT prob SE df lower.CL upper.CL
## none 0.988 0.0209 79 0.7055 1.000
## crabs 0.751 0.3021 79 0.1081 0.987
## shrimp 0.614 0.3790 79 0.0619 0.975
## both 0.325 0.3458 79 0.0204 0.917
##
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale
##
## $contrasts
## contrast odds.ratio SE df t.ratio p.value
## SYMBIONT.crab_vs_shrimp 1.89 2.18 79 0.556 0.5800
## SYMBIONT.one_vs_both 4.55 4.75 79 1.454 0.1499
## SYMBIONT.symbiont 62.91 79.54 79 3.276 0.0016
##
## Tests are performed on the log odds ratio scale
emmeans(mckeon.glmmTMB1, ~SYMBIONT) %>% regrid() %>% contrast(list(SYMBIONT=cmat))
## contrast estimate SE df t.ratio p.value
## SYMBIONT.crab_vs_shrimp 0.137 0.254 79 0.539 0.5916
## SYMBIONT.one_vs_both 0.358 0.225 79 1.589 0.1161
## SYMBIONT.symbiont 0.425 0.296 79 1.437 0.1545
Conclusions:
r.squaredGLMM(mckeon.glmmTMB1)
## Warning: The null model is correct only if all variables used by the original
## model remain unchanged.
## R2m R2c
## theoretical 0.1489321 0.8674549
## delta 0.1437647 0.8373578
## The delta method can be used with for all distributions and link
## functions, while lognormal approximation and trigamma function are
## limited to distributions with logarithmic link. Trigamma-estimate
## is recommended whenever available. Additionally, for binomial
## distributions, theoretical variances exist specific for each link
## function distribution.
performance::r2_nakagawa(mckeon.glmmTMB1)
## # R2 for Mixed Models
##
## Conditional R2: 0.867
## Marginal R2: 0.149
emmeans(mckeon.glmer1, ~SYMBIONT, type='response') %>%
as.data.frame %>%
ggplot(aes(y=prob, x=SYMBIONT)) +
geom_pointrange(aes(ymin=asymp.LCL, ymax=asymp.UCL))
emmeans(mckeon.glmmTMB1, ~SYMBIONT, type='response') %>%
as.data.frame %>%
ggplot(aes(y=prob, x=SYMBIONT)) +
geom_pointrange(aes(ymin=lower.CL, ymax=upper.CL))
# References